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We perform two flavor QCD simulations with an imaginary chemical potential and measure 
derivatives of the pressure up to 4th order as a function of the imaginary chemical potential 
and the temperature T G [0.83T C ,2T C ]. For temperatures T > T c , these derivatives are fitted by 
a Taylor series in fi/T about ji = 0. A fit limited to 4th order describes the data poorly at all 
temperatures, showing that we are sensitive to 6th order contributions. Similarly, a 6th order fit 
fails for temperatures T c <T < 1.05X C , showing the need for 8th order terms. Thus, our method 
may offer a computational advantage over the direct measurement of Taylor coefficients at jj. = 0. 
At temperatures T < T c , we fit our data with a hadron resonance gas ansatz. The fit starts to fail 
at T > 0.9571-. Using our fits, we also reconstruct the equation of state as a function of real quark 
and isospin chemical potentials. 
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1. Introduction 

Although lattice QCD has been used successfully for simulations at zero and finite temper- 
atures and at zero density, Monte Carlo simulations at non-zero densities suffer from a technical 
problem: the lattice QCD action becomes complex, which prevents its customary probabilistic 
interpretation. In principle one could perform simulations at zero density, and use the reweight- 
ing technique to obtain information at finite densities. An early attempt known as the Glasgow 
method [jl|] did not work due to the overlap problem: the configurations at zero density were too 
"far" from the target configurations at non-zero densities. Considerable progress has been accom- 
plished by generalizing the Glasgow method to two-parameter reweighting Nevertheless, the 
range of reliability of this technique is difficult to assess, and its failure can go undetected. 

Therefore, another, more conservative way to deal with finite baryon densities may be useful. 
It consists of calculating Taylor coefficients of observables with respect to the chemical potential 
/i about jU = 0. Those Taylor coefficients can be expressed as expectation values of complicated 
observables, which can be measured at zero density. Thus, there is no difficulty to perform Monte 
Carlo simulations in this method. A first, pioneering attempt to obtain quark susceptibilities [||] 
has been followed by numerous works, obtaining in particular the response of screening masses to 
chemical potential [|], |5| ^J. The Taylor expansion method has also been used for studies of the 
equation of state, of the phase transition and of higher order susceptibilities [Q, [g, ||, 10]. However, 



the complexity of the observable representing the Taylor coefficient, and the computer effort to 
measure it, increase rapidly with the order of the Taylor expansion. This motivates us to follow a 
different strategy. 

Since no difficulty appears for simulations at imaginary chemical potential jj, = ifa, one can 
obtain information at finite baryon densities by analytic continuation of observables measured at 
finite fa. Actually, this imaginary chemical potential strategy has been applied with success to the 
determination of the phase transition 

In this study, we perform simulations at finite fa and measure derivatives of the pressure as 
a function of fa. These derivatives contain information about the Taylor coefficients of the /I = 
expansion, which can be extracted by fitting. Finally, we try to reconstruct the equation of state at 
finite baryon and isospin densities. The strategy of our method and preliminary results were pre- 



sented rather long ago in [ 12]. Here we report further progress on this project. A related approach, 



where the quark density is measured at imaginary quark and real isospin chemical potentials and 



then fitted by a polynomial ansatz, has recently been presented in JT5[|. 

2. Equation of State at Finite Chemical Potential 

The lattice QCD partition function with Nf flavors of staggered fermions can be written as 

Z = Jnf f detM(U,m h Hi) 1/4 ex V (-S g [U])dU, (2.1) 

where S g [U] is the gauge action and M(U, nit, fa) stands for the staggered Dirac operator with quark 
mass m; and chemical potential fa . In this study we consider Nf = 2 degenerate fermion species 
and use the standard Wilson gauge action. 

'We set aside potential problems with "rooting" the determinant, particularly at non-zero chemical potential. 
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The pressure or the equation of state with chemical potential \i u and \l d is given by 

F T 

p(Li u ,Hd) = -y = yfoZ(Ll u ,Hd), (2-2) 
and can be expanded in a Taylor series about \i u = ji^ = as 

Ap _ p(^,ii d )-p(0,0) _ ^ 1 (t*Y(t£\ n (2 3) 



T 4 T 4 , n\m\ \T J \T 

n,m=l 

where f nm are the Taylor expansion coefficients. They vanish when (n + m) is odd due to CP 
symmetry. Furthermore, for equal quark masses there is another symmetry f nm = f mn . The f nm 's 
are related to derivatives Xij of the pressure measured at non-zero chemical potential by 

T i+r-A Y _ d i+ J(p(n u ,H d )/T*) _ y 1 (toy-* (24) 

XlJ d{pL u /TYd^ d /T)i n= t =j {n-i)\(m-j)\ }nm \T) \T ) ' ^ 

While at zero density Xij = fijT 4 ~'~ J , at non-zero densities Xij includes higher order /,„„ terms, 
and does not vanish for odd (i + j). This suggests to use all available Xij's at non-zero densities, 
in order to estimate the f nm 's. Here, we try to estimate by fitting all %ij simultaneously to the 



polynomial expansions eq.(|2.4]). Of course, Xij at non-zero baryon density is not directly obtainable 
from simulations on the lattice because of the sign problem. However, Xij can be obtained through 
simulations at imaginary quark chemical potential or at real isospin density. Here, we calculate %ij 
at imaginary chemical potentials. 

Therefore, we set ji. = ijij. Each Xij depends on higher order Taylor coefficients following 



eq.([L4|). Therefore, with sufficiently accurate data on %ij one can a l so obtain higher order Taylor 
coefficients f nm , n > i,m > j. The measurements of the derivatives involve computing traces of 
inverse Dirac matrix products. These traces were estimated using the noise method with 40 Z2 
random vectors. In this study we measure Xij U P to i + j = 4. Thus we have 8 different Xi/ S - We 



fit all the data to the corresponding 8 polynomial expansions eq.(p.4|) truncated to a given order 
[n + m), and try to obtain the Taylor coefficients f nm . 

As we will see, in the confined phase a Taylor expansion is not the most compact description 
of the pressure. Instead, for T < T c we use the Hadron Resonance Gas (HRG) model. In the HRG 
model the pressure is given as 2 

= G[cosh(^)-l] + / ? [cosh(^)cosh(f )-l] (2.5) 
+ W[cosh(^i) (cosh(^)+cosh(^)) -2], 

where G,R and W are constants related to the hadron spectrum, and quark and isospin chemical 
potentials ji q and jij s are defined as \i q = + Hd)/2 and jij s = (jl u — Hd)/2 respectively. 

The derivatives of the pressure with respect to \i u and \l d , instead of having the polynomial 



form eq.(2.4), are now obtained by differentiating eq.(2.5). The coefficients G,R and W are then 



"This expression is taken from (4.3) in [|8h 
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Figure 1: % 2 jdof of various polynomial ansatze as a function of T/T c . The fitting range of ajij is 0.0-0.24. 



extracted by fitting imaginary-/^ data. In terms of G,R and W, the first Taylor coefficients are given 

(2.6) 

(2.7) 
(2.8) 

(2.9) 
(2.10) 



f20=[G + -R + lwy 

fn = -(G + 2(R+W)), 
f 22 = G + 4(R + W), 

f il =G + 1 -R + 49W, 
f 40 = -G + 5(R + W). 



Table 1: % jdof for polynomial ansatz of degree 4, 6 and 8 (maximum value of (n + m) in eq.(2.4)) 



T/T c 


0.99 


1.00 


1.03 


1.04 


1.065 1.085 


1.1 


1.2 


1.3 


1.4 


1.5 


2.0 


4th 


85.1 


134.9 


3.15 


3.14 


3.50 7.24 


3.20 


10.9 


11.3 


6.67 


9.95 


9.15 


6th 


19.1 


42.1 


1.60 


2.19 


0.82 5.50 


5.53 


1.52 


0.89 


2.46 


1.09 


2.10 


8th 


4.53 


5.29 


1.64 


1.77 


0.81 1.01 


2.15 


1.72 


0.91 


2.16 


1.22 


1.28 



Table 2: % 2 /dof for HRG ansatz. 



T/T c 


0.83 


0.9 


0.95 


0.98 


0.99 


1.00 


HRG 


1.29 


1.00 


2.10 


15.8 


10.5 


29.4 



3. Simulations at Imaginary Chemical Potential 

We have performed simulations on 8 3 x 4 lattices at a quark mass m q = 0.05 and imaginary 
chemical potentials a\i\ = 0.0, . . . ,0.24. We have chosen 16 values of /3 ranging from 4.90 to 6.85, 
which correspond to T /T c = 0.83 ~ 2.0. Most of the simulations were performed using the R- 
algorithm with a step size At = 0.02. We also used the Rational HMC algorithm [13] to check 
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(a) (b) 

Figure 3: Equation of state (Pressure) as a function of a\i q and a\ij s at (a) j3 = 5.53 (T /T c ~ 1.1) and at 
(b)j3 = 4.90(7/% ~ 0.83). 



the systematic stepsize errors caused by the R-algorithm, and found no significant difference for 
this lattice size and quark mass. At each simulation point we have accumulated 12000 to 20000 
measurements. The measurements were taken every 5 trajectories to balance the computational 
effort of the R-algorithm simultion and measurements. 



3.1 Fitting to Xij 

We determine f nm by fitting all the derivatives simultaneously to the corresponding ansatz of 



Xij. We used the polynomial ansatz eq.{2A) for the data at T/T c > 0.99, and the HRG ansatz 
eqX^) at T /T c < 1.0. Tables 1 and 2 show the % 2 jdof for the polynomial and HRG fits, respec- 
tively. The fitting range of a\i\ is 0.0 — 0.24, which covers most of the range up to the Roberge- 
Weiss transition at \i\ = 7iT /3. 

Fig.|l| compares the x 2 /dof among various polynomial and HRG fits. One can see that the 4th 



order polynomial ((n + m) < 4 in the expansion eq.(2.4)) is not good over the whole temperature 
range, and that the 6th order one becomes poor in the vicinity of T c . Similarly, one can also see 
that the quality of the fit based on the HRG ansatz becomes poor for T /T c > 0.95. While the 



failure of the HRG ansatz near T c has been noticed before [14, 15], it is remarkable that we can see 
clear indications of 6th order, and even 8th order Taylor coefficients with our modest study. The 
measurement of 8th order Taylor coefficients represents the current state of the art JIo|]. 

Fig.|| shows the Taylor coefficients fu, fzo, fn, fa and fa as a function of temperature. 
Those results are obtained by fitting a 6th order polynomial in a range of a/j,j = 0.0 — 0.24. For 
T/T c > 1.0, they agree well with those obtained from the direct measurement of derivatives at 
/I = 0, i.e. Xij\u=0i b ut 316 more accurate. We do not show the 6th order Taylor coefficients 
fa: fa -,3 42 and ^33: even though their collective effect is statistically significant, they cannot be 
individually determined with any statistical accuracy. We only observe that fa is dominant at this 
order. 

Similary, the Taylor coefficients obtained from the HRG ansatz for T /T c < 0.95 also agree well 
with direct measurements of Xij\u=o> w i tn higher accuracy. However, for T /T c > 0.95 the results 
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(a) (b) 
Figure 4: Number density at j3 =5.53(T/T C — 1.1) as a function of a\l q anda^/ s : (a)N q /T 3 and (b)Ni s /T 3 . 




(a) (b) 

Figure 5: Number density at j8 = 4.90(T/T C ~ 0.83) as a function of a\i q and a\i h : (a)N q /T 3 and (b) 
N Is /T\ 



from the HRG ansatz fits deviate from Xij\fi=0- This observation is consistent with the measured 
X 2 /dof, which increase considerably for T/T c > 0.95. 



3.2 Equation of State at Finite Densities 

Once we obtain the Taylor coefficients of the pressure or the parameters of the HRG model, 
we can reconstruct the equation of state. Here, we present two cases at j8 = 5. 53 (T/T c ~ 1.1) 
and j3 = 4.90(T /T c ~ 0.83) which are reconstructed with the Taylor series and the HRG ansatz, 
respectively. Fig. ||(a) shows the equation of state at j8 = 5.53(T/r f ~ 1.1) as a function of a\i q 
and a\Xi s . Similarly Fig. |(b) shows the equation of state at j3 = 4.90(r /T c ~ 0.83). 

One can also reconstruct other interesting quantities. Fig. ^ shows the quark number density 
N q and the isospin number density Nj s at j8 = 5.53 as a function of a\l q and a/j.j s . Similarly, Fig. || 

shows N q /T 3 and N Is /T 3 at j3 = 4.90. Here, N q and N Is are defined as N q = and N Is = 

' ' o\l q d\li s 

respectively. 
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4. Conclusions 

We have performed simulations at imaginary chemical potentials and measured the derivatives 
of the pressure with respect to pL, at zero and non-zero imaginary /x. By fitting all the derivatives to 
a polynomial ansatz or an HRG ansatz, we obtained the Taylor coefficients of the \i/T expansion 
of the pressure about /i = 0. The Taylor coefficients obtained by a polynomial fit for T/T c > 1.0 
agree well with the direct measurement of derivatives at /i = 0, Xij\ix=Q-> but are more accurate. 
Remarkably, we find it impossible to obtain a good fit, at any temperature, without including 6th 
order derivatives. For T c <T < 1.04r o 8th order derivatives are necessary. Thus, our approach 
may provide a cheaper alternative to the direct measurement of high-order derivatives at fL = 0. 

Similarly, below T c we observed that the Taylor coefficients obtained by the HRG ansatz de- 
viate from Xij\u=o f° r T /T c > 0.95, and the HRG ansatz itself gives a poor description of the 
imaginary-jU data. The same observation has been made in [|l5|]. 

Finally, using the obtained Taylor coefficients we reconstructed the equation of state and the 
number densities as a function of \i q and \i\ s up to 4th order. 
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